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Abstract: We present results of numerical simulations of the Nambu - Jona- 

Lasinio model with a non-zero baryon chemical potential //, with particular emphasis 
on the superfluid diquark condensate and associated susceptibilities. The results, when 
extrapolated to the zero diquark source limit, are consistent with the existence of a 
non-zero BCS condensate at high baryon density. The nature of the infinite volume and 
zero temperature limits are discussed. 
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1 Introduction 

o 

Colour superconductivity (CSC) at low temperature and high baryon number density is 
one of the most intensively studied topics in QCD thermodynamics (for recent reviews 
see IJ). The basic idea is that diquark pairs at the Fermi surface in quark matter be- 
come bound and condense in a relativistic analogue of the BCS mechanism for electronic 
\ superconductivity @. Because of the strong attractive qq force in QCD, however, es- 
timates of the energy gap A, which opens up at the Fermi surface, can be as large as 
lOOMeV ||, with potentially important consequences for the physics of compact stars. 
For instance, one intriguing possibility is that pairing between quarks with differing 
Fermi momenta may lead to a crystalline phase in which translational invariance of the 
ground state is spontaneously broken [||. 

Analytic approaches to the problem, however, are restricted either to asymptotically 
high (and phenomenologically irrelevant) densities where perturbative QCD is applica- 
ble, or to self-consistent treatments of model theories which capture some but not all 
relevant physics. A first principles calculation using lattice QCD remains elusive due to 
the difficulties of performing Monte Carlo simulations with baryon chemical potential 
/i ^ 0. There are, however, simpler models which are simulable. Two Colour QCD is 
a confining theory in which the lightest baryons are bosonic qq states. At high density 
its ground state has a gauge invariant condensate (qq) ^ which spontaneously breaks 
baryon number symmetry leading to superfluidity. Because of the absence of a Fermi 
surface, however, the phenomenon far more closely resembles a standard Bose-Einstein 
condensation, and is well-described by chiral perturbation theory ||; this picture is 
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confirmed by lattice simulation [H, M. The NJL model M, by contrast, is a strongly- 
interacting model with fermionic excitations but no confinement. It can be treated by 
self-consistent analytic methods and has been applied to CSC || as well as more general 
issues in low-energy QCD ||. 

The NJL model with /i 7^ has also been simulated, on a 2+1 dimensional lattice. 



Although evidence for qq pairing has been found in the scalar isoscalar channel [10 



BCS condensation does not take place ||11|| , the argument being that long- wavelength 



fluctuations in the phase of the condensate wavefunction wipe out long range order. 
Because phase coherence remains, however, superfluidity is realised in an unconventional 
"thin film" fashion. This picture, related to the low dimensionality of the system, is 
essentially non-perturbative and would not be exposed by use of self-consistent methods. 
It points to the potential importance of fully non-perturbative treatments of even simple 
models in attempts to understand CSC. In this Letter we extend the lattice analysis of 
[IT|| to the NJL model in 3+ld, motivated by (i) the possibility for a genuine BCS 



condensation in this case, and (ii) the model's phenomenological relevance. In the 
following sections we review the lattice formulation of NJL, show how lattice parameters 
may be chosen to match low energy QCD, and then report results from simulations 
with fi 7^ 0, paying particular attention to diquark observables. We will show that 
the high density phase is qualitatively different to that observed in 2+ld and far more 
closely resembles a conventional BCS superfluid; however, some issues concerning the 
thermodynamic limit need to be resolved before this conclusion can become definitive. 



2 Lattice Model and Parameter Choice 



The model studied here is the 3+ld version of that studied in JTTJ. In particular it has 
the action 
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S = V tr AV + — E l^ 2 + ™ 
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(1) 



where g 2 is the four-fermi coupling constant, the bispinor ^ is written in terms of 
staggered isospinor fermionic fields via ^f tr = (%, x* r ), and the auxiliary bosonic fields 
a and 7?, which live on sites a; of a dual lattice, are introduced in the standard way. 
Written in the Gor'kov basis the fermion matrix is 



A 



1 

2 



jr 2 M 
-M tr JT 2 



(2) 



where the matrix M is defined by 
1 
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( eM( Wo - e ^yx-b) + {5 yx+0 - 5 yx „ p ) + 2m 5 xy 



(x,x) 



(3) 
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mo is the bare quark mass and the symbols i] u (x), e(x) are the phases (— i) x o+---+^-i anc [ 
(_]_ ^0+^1+^2+2:3 respectively. The Pauli matrices f are normalised such that tr(r i r 7 ) = 
25ij, and (ic, x) represents the set of 16 dual lattice sites neighbouring x. The diquark 
sources j and j are introduced to allow us to measure the diquark condensate and play 
the same role as a bare quark mass in the lattice estimation of the chiral condensate. 
(NB. The sources used are greater than those in [JlTJ by a factor of 2, which allows 
us to identify them with Majorana fermion masses.) We simulate the action ([!]) using 
a hybrid Monte Carlo (HMC) algorithm which uses a functional weight det^4^4^, thus 
doubling the number of fermion degrees of freedom, but which has the advantage of of 
being exact. We treat the diquark terms in the partially quenched approximation, ie. 
the sources j and j are set to zero during the HMC update of the bosonic fields, but are 
made non-zero during the measurement of diquark observables. Field theoretically, this 
means that the formation of virtual diquark pairs in the vacuum is suppressed. This 
approach has the advantage that one can examine many source strengths for only one 
(computationally expensive) chain of field evolutions. 

Our formulation employs staggered quarks, which, although preserving some of the 
model's continuum chiral symmetry, increase the number of quark degrees of freedom by 
a factor N c = 4. As the NJL model has no gauge degrees of freedom we interpret these 
doublers phenomenologically as "colours". The model is therefore an effective Nf = 2, 
N c = 4 QCD-like theory. It should be stressed, however, that a diquark condensate in 
this model breaks global rather than local symmetries, hence physically giving rise to 
superfluid rather than superconducting behaviour. The continuum limit of the model 
yields the Lagrangian density |T2], Q 

- (Vi75 ® T ® 75^i 

1 r 
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£cont = Ipi W + m o + ^70) A ~ y 



+ 2 



j (i>i tr Cj 5 <g> r 2 <g> C 7 5^) + J (V^7s ® r 2 <g> C 75 ^/ r )] , (4) 



which in the limits m — > and j, J — > recovers the SU(A^j)^ ® SU(Nf)ji chiral and 
U(1)b baryon number symmetries of QCD respectively^. In the tensor products the first 
matrix acts on spinor, the second on isopinor, and the third on "colour" indices. 

Due to the point-like nature of the four-fermi interaction, the NJL model has no 
interacting continuum limit for d > 4, leaving physics dependent on the UV regularisa- 
tion employed 0, l^j. In order to ensure we simulate in a physically plausible regime, 
we employ the methods used in || to match our model's parameters to low energy, 
vacuum phenomenology, taking advantage of the fact that a perturbative expansion in 
1/N C is possible in four-fermi theories. We calculate quantities analytically to leading 
order in 1/N C , the Hartree approximation, using staggered quark propagators defined 
on a L? s x L t /2 4 Euclidean blocked lattice with periodic boundary conditions in 
spatial dimensions and anti-periodic boundary conditions in the temporal dimension. 
The momentum loop integrals are discretised using mode sums and taken to the limit 
V- 1 -> where V = L 3 L t . 



1 Doubling due to det^4^.T results in an additional SU(2AT C ) global symmetry in the continuum limit. 
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In particular we calculate the dimensionless ratio between the pion decay rate f n and 
the constituent quark mass m* = £ + m , where £ = (a). 



f d 4 p 



cos 2^ 



f (2tt) 



p 2 + (am*) 2 \ 



1 



-I (2tt) 



p 2 + (am*)' 



2 ■ 



(5) 



where p 2 = X^=o s i n2 p v . 

Fixing f n to its experimental value of 93MeV and m* to a physically reasonable 
350MeV, we extract a dimensionless quark mass of am* = 0.3251, where a is the spacing 
between lattice points, which in this section we will not set to unity. Calculating the 
product of £ and the dimensionless inverse coupling (3 = a 2 /g 2 , 



Epa = 2N c N f am* J\ 



I d 4 p 



(27r) 4 p 2 + (am*) 2 



0.1826, 



(6) 



allows us to derive that (3 = 0.1826/(0.3251 — amo). 

Finally, calculating the mass of the pion m n allows us to fix mo, and therefore (3. In 
the Hartree approximation the pion propagator can be written as a series of connected 
bubble diagrams 



g 2 g 2 Tl ps (k 2 )g 2 g 2 Il ps (k 2 )g 2 Il ps (k 2 )g 2 
k 2 + (am n ) 2 2 4 8 

in terms of the vacuum polarisation of the pion 
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2 + 2 2 n ps (P) 



(7) 



n ps (P) 



m — mo 

m*g 2 



+ 2N c N f k -l(k 2 ), 



where 
l(k 2 



d 4 p 



f (2tt) 



J2fj, sin 2 (p + |) + (am*) 2 J2„ sin 2 (p - |) + (am*) 



(8) 



(9) 



p^ is the loop momentum shifted to be independent of a, and k^/a is the physical 
momentum of the pion. The fact that there is a pole in (|7]) at k 2 = —a 2 m 2 combined 
with (0) allows us to write 
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% 

™* N c N f g 2 l(~k 2 



(10) 



A bare quark mass of amo = 0.002 leads to a phenomenologically acceptable pion 
mass of 138MeV and implies an inverse coupling (3 = 0.565. Finally, as we know both 
am* and m* we can extract the lattice spacing in the large- N c approximation, with the 



result a 



1076MeV, or a ~ 0.2fm. 
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3 Equation of State 



As has been previously discussed, in the limit mo — > the NJL model has an exact 
SU(2)l®SU(2)/j chiral symmetry, which for sufficiently large g 2 is spontaneously broken 
in the vacuum to SU(2) isospin , leading to a dynamically generated quark mass m* = S, 
and 3 massless Goldstone pions. In the presence of a baryon chemical potential /x the 
symmetry is approximately restored as /x is increased through some onset scale /x Q ~ S, 
with the order of the transition being highly sensitive to the parameters employed ||. 

We determine the nature of this transition in our physically reasonable regime by 
studying the order parameter, the chiral condensate (xx)i defined by 

1 dlnZ . . 

<»> = vS^- < n > 

where Z is the partition function. We also study the baryon number density n#; 

nB = w^r- (12) 

These are calculated as functions of /x, both in the large-iV c limit and by numerical 
lattice simulation, using a repeated stochastic estimator for the diagonal elements of the 
inverse fermion matrix. We use a HMC algorithm with (3 = 0.565, mo = 0.002 on x L t 
lattices with L s = L t = 12, 16 and 20 for various values of /x, as well as on lattices with 
L s L t at fx = 1.0. Approximately 500 equilibrated configurations separated by HMC 
trajectories of mean length 1.0 were generated in each run, with measurements carried 
out on every other configuration. Statistical errors were calculated using a jackknife 
estimate. 

Our results are shown in Fig. [l|, where the solid curve denotes the solution in the 
Hartree approximation, where (xx) = and the points denote lattice results. We 

carried out a linear extrapolation in V" 1 to the thermodynamic limit. To leading order 
in 1/N C chiral symmetry is approximately restored via a crossover between 0.4^ii^0.6. 
The lattice data agree qualitatively with this although both (xx) an d fJ> are about 30% 
smaller, which we attribute to 0(1/N C ) corrections, ub increases approximately as a 
power of fi. 

4 Diquark Condensation 
4.1 Observables 

In order to explore the possibility of a U(1)b- violating BCS phase at high /x we study 



both the diquark order parameter and susceptibilities |TT|. The operators 



qq±(x) = x tr jx(x) ± xjx tr (x) (13) 
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Figure 1: Chiral condensate and number density as functions of /i extrapolated to V' 
showing both the large- N c solution (solid curve) and lattice results (points). The 
diquark condensate is plotted in the zero temperature, zero-source limit. 



allow us to define the diquark condensate as 



-) 



1 d\nZ 



V dU ' 

where j± — j ± j. We also define the susceptibilities 

1 



(14) 
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E {X tr T2X(0)X tr T2X(x) +XT2X tr (0)XT2X tr (x] 

x 

± (x tT T2X(0)XT2X tr (x) + XT2X tr (0) X tr T2X(x)) . (15) 



These can be expressed as the sum of two connected contributions corresponding to the 
two possible Wick contractions 



X 



(trrg xx ) 2 ) - (trreu 2 ] + (trg 0x rg£) = x s + x T: 



(16) 



where Q = A^ 1 is the Gor'kov propagator and T projects out the appropriate compo- 
nents. By analogy with meson physics we label these components "singlet" and "non- 
singlet" respectively. Both (qq + ) and x± are calculable using the same stochastic esti- 
mation method used to measure (xx) an d n B , whereas x± s is measured using standard 
lattice methods for meson correlators. We fix j = j = j* throughout. It is interesting 
to note that although in most cases the singlet contributions are found to be consistent 
with zero, in the low /i phase with large j, x+ can be up to 10 — 20% the magnitude 
of x+ s - Therefore in contrast with the NJL model in 2+ld [^T|] we cannot ignore these 
contributions and assume that x+ — XT- 
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From ( P^D it is straightforward to derive the Ward identity 



which along with the ratio 



R = -%±- (18) 

allows one to distinguish between the two phases as j — > 0. With U(l)# manifest, the 
two susceptibilities should be identical up to a sign factor and R should equal 1 in this 
limit. If the symmetry is broken, the Ward identity predicts that \- should diverge and 
R vanish. 



4.2 Results 



The susceptibilities were measured and R calculated for the aforementioned lattice sizes 
at various values of fi. We extrapolated the data using least-square linear fits to the 
limit L^ 1 — > 0, which experience in 2+ld |TT] has shown to be the best model of finite 



size corrections. 
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Figure 2: i? as a function of j for various ji. 



Fig. H shows this extrapolated data plotted against j for various values of \i. One can 
immediately notice that although a linear fit through the data for j > 0.3 is reasonable, 
for j < 0.3 the data departs sharply from the fit, especially in the high density phase 
/x~0.6. Possible origins of this effect are explored in section fO|. Assuming that we can 
disregard the points with j < 0.3, one can see from Fig. |2] that for /i = a linear fit is 
consistent with a ratio of R ~ 1 corresponding to a manifest baryon number symmetry 
as one would expect in the vacuum. As /i is increased the value of the ratio decreases 
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and as fi approaches one inverse lattice spacing we see that R ~ 0, suggesting that the 
U(1)b symmetry is broken. 

For more direct evidence of diquark condensation we study the diquark condensate 
defined in ([14]). Fig. ^| shows (qq+) plotted against j for various values of /i. Again we 
have extrapolated linearly to the limit L^ 1 — > 0. We use a quadratic fit through the 




0.2 0.4 0.6 0.8 1 

j 



Figure 3: (qq + ) as a function of j for various fi. 

data with j > 0.3 but again one can clearly see that for high /i the low j points lie well 
below the curves. Ignoring these points we extrapolate to j — > 0. For = we find 
no diquark condensation as expected, but as /i increases from zero, so does (qq+). We 
believe that the observations lim^o-R = and lim^o 7^ together support the 
existence of a BCS phase at high chemical potential. 

Finally, (qq+) is plotted as a function of fi in Fig. [I]. Although there is clearly 
a crossover from a phase with no diquark condensation to one in which the diquark 
condensate has a magnitude approximately that of the vacuum chiral condensate, this 
crossover is far less pronounced than in the chiral case. (qq+) increases approximately 
as /j, 2 , (ie. as the surface area of the Fermi-sphere), but eventually saturates and even 
decreases as /i increases through 1. This may be a cut-off artifact. The curvature 
d 2 (qq+) /dfi 2 is positive, in contrast to the behaviour observed in simulations of Two 
Colour QCD in which there is no Fermi surface and U(l)s breaking proceeds via Bose- 
Einstein condensation 0. It is possible, of course, that this behaviour for intermediate 
fi is an artifact of our poor control over the j '• — > extrapolation. 

4.3 The Low Temperature Lattice Fermi Surface 

We know, from experience in 2+ Id, that an extrapolation in L^ 1 is the best description 
of of finite size corrections in this model. One natural progression would be to simulate 



S 



with a small L s and various L t , saving both CPU time and resources. However, this 
approach presents a problem. 

As we increase \x we build up a Fermi surface, which in the continuum is a sphere 
smeared out by thermal fluctuations over a region 5k F ~ 2T = (2L t )~\ In Euclidean 
space, an increase in L t corresponds directly to a decrease in temperature, and at low 
temperatures, the Fermi-Dirac distribution closely resembles a step-function with all 
states with sinh -1 \Jj2l=i sin 2 k < Ep occupied, and all other states unoccupied. If we 
simulate with L t 3> L s , the smearing of the Fermi surface will be too fine to be resolved 
on the coarse momentum lattice, and as \x ~ Ep is increased the physics will be constant 
except when the Fermi sphere crosses each momentum mode, turning any transition into 
a series of steps. 




Figure 4: Chiral condensate and number density in the large- N c limit for various lattices. 

Fig. |] illustrates the effect that this has on the chiral transition in the large- N c limit. 
(xx) an d i^b are plotted on 12 4 , 12 3 x 20 and 12 3 x 40 lattices. On the asymmetric lattices 
one can clearly see the discontinuities caused as the Fermi surface crosses momentum 
modes. In the baryon number density one can even observe a plateau as the first mode 
becomes occupied. This behaviour was also observed in some initial lattice simulations. 
In practice, it seems that the effect is only prominent when L t ^3L s /2. 

4.4 Finite Volume Effects 

To investigate the possibility that the low-j discrepancies in Figs. || & [| are due to finite 
volume errors, it is necessary to study spatial volume dependence in a controlled manner. 
We study all L 3 x L t lattices with L s = 6, 8, 10, 12, 16 & 20 and L t = 12, 16 & 20 
for fixed \i — 1. Being both well into the high-/i phase and far from the transition, any 
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finite volume effects should be easily identifiable. We extrapolate to infinite temporal 
extent, corresponding physically to zero temperature, using the fact that the dominant 
finite size correction in our data is proportional to L^ 1 . However, Fig. [5] illustrates that 
there is some residual L s dependence. 





o 




















Figure 5: (qq+) in the zero-temperature limit as a function of the spatial extent of the 
lattice. The results are scaled such that {qq+)\j 



=20 



1. 



(qq+) is plotted as a function of spatial extent for various values of j, with data 
presented as fractions of (qq+)\ L =20 to enable them to be displayed on the same axes. 
For j > 0.5, L s > 12 the effect is negligible; the data are consistent within errors. 
For j < 0.5 we observe non-monotonic fluctuations, whose magnitude increases with 
decreasing j. A possible explanation is that with a small source, the region about the 
Fermi surface in which diquark pairs would bind may be too fine to be resolved on 
small spatial lattices. Variational studies of the Nf = 2, N c = 3 continuum NJL model 
at zero temperature in a finite spatial volume |T3| find that with no diquark source, a 
spatial extent of 7fm (~ 35 lattice spacings) is required before the model approximates 
its infinite volume limit. Below this volume, the BCS gap oscillates rapidly about its 
thermodynamic limit solution. Although this may explain why, with a small source, the 
condensate displays this non-monotonic L s dependence, for L s > 12 these fluctuations 
are less than a 1% effect, whereas the depletion of diquark pairs below the extrapolated 
curve in Fig. |3] is ~ 30%. At this stage, therefore, we cannot claim to understand the 
origin of the depletion of diquarks at j < 0.3. 



5 Summary and Outlook 

We have presented evidence for a BCS phase in the 3 + Id lattice NJL model at high 
chemical potential, by measuring a diquark order parameter which appears, at n — 
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1, to be of similar magnitude to the chiral order parameter in the vacuum. This is 
supported by a susceptibility ratio, lirn^o R(fi = 1) « 0, which suggests a broken U(l)# 
baryon number symmetry in the high-/z phase. Both features contrast sharply with the 
equivalent measurements in the 2+ld NJL model where (qq+) vanishes as a power 
of j and R is approximately a constant independent of j, both indicative of a critical 
phase with unbroken U{1)b symmetry at high density. 

Our measurement of (qq+) can be translated into an order-of-magnitude estimate for 
the gap A using the following simple-minded argument: the condensate is the number 
density of bound diquark pairs, which in a BCS scenario is roughly equal to the volume 
of fc-space within a distance A of the Fermi surface divided by the density of states; this 
yields (qq + ) ~ A/i 2 . Since at /ia = 1.0 we find (qq+)a 3 ~ m*a ~ 0.2 — 0.3, it is quite 
plausible that the gap values of O(100MeV) predicted in || are possible. 

Our claim to have observed a BCS phase in a lattice simulation is necessarily pro- 
visional, however, because it depends on the discarding of data with j < 0.3, which in 
the high /i phase fall rapidly away from the trends observed at higher j. This threshold 
is very high when compared with chiral symmetry breaking at « = 0, where the bare 
quark mass can be taken as low as m = 0.002 with no adverse effects on (xx) on the 
volumes studied. It is also larger by a factor of 3 than the equivalent threshold observed 



in 2+ld [11]. The volume effects discussed in section WA do not resemble those expected 



from the standard treatment of Goldstone fluctuations [15|]. Variational studies in fl4j 



suggest that with a vanishing source, a 7fm spatial box (L s ~ 35a), which is currently 
unattainable, may be required for the BCS state to approximate its infinite volume limit. 

It is also possible that the low-j discrepancy is a result of the partially quenched 
approximation, and would be resolved in a full unitary simulation. Although a repeat of 
the current study with a full simulation is currently beyond our reach, were it deemed 
necessary, it would be sufficient to repeat the study presented in section JO] with dynam- 
ical diquarks for a limited range of j. It should be noted, however, that no significant 
effect on diquark observables due to partial quenching was noted in flTT|. 

We have chosen a value m* = 350MeV for the constituent quark mass in an attempt 
to simulate with phenomenologically reasonable parameters. The study of was done 
with m* set to an equally reasonable 400MeV. Initial studies of the large-iV c lattice 
model have determined that with this choice, the chiral transition changes character to 
become much sharper and may even become first order in the infinite volume limit. It 
is an interesting coincidence that this change in the order of the transition may occur 
in the "physical" region of parameter space. Not only is this worth exploring, but it 
might make the determination of a BCS phase more clear cut. It is also desirable to 
determine the lattice spacing a by direct measurement of fermion and pion masses so that 
we can convert our results into physical units independent of the large- N c assumption. 
Finally, we wish to study the fermion dispersion relation, which will allow us to compare 
directly the chiral mass gap S with the BCS gap A, these being, in principle, physically 
measurable quantities. 
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